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INTRODUCTION 


The  ratios  of  spherical  Bessel  functions  of  complex  argument  are  needed  for 
Mie  calculations  of  scattering  and  absorption.  Lentz  has  developed  a  new 
continued  fraction  technique1  for  evaluating  these  ratios,  with  an  improvement 
to  eliminate  errors  when  certain  terms  approach  or  become  zero.  The  authors 
of  a  recent  note,1  however,  had  difficulty  in  applying  the  algorithm 
improvement  when  either  the  numerator  or  denominator  was  exactly  zero.  This 
case  is  special,  and  a  simplification  is  possible.  The  purpose  of  this  paper 
is  to  present  an  algorithm  for  this  special  case,  together  with  some 
clarifying  remarks  concerning  its  development. 


BACKGROUND 

The  Lentz  algorithm  is  simply  a  way  of  calculating  continued  fractions;  the 
algorithm  begins  at  the  beginning  of  the  continued  fraction  rather  than  at 
some  indeterminate  final  coefficient.  An  improvement  formula  allows  one  to 
bypass  any  step  in  which  a  given  numerator  or  denominator  might  be  zero  (to 
the  accuracy  of  the  computer).  In  addition,  roundoff  errors  may  be  reduced  by 
using  the  algorithm  improvement.  When  applied  to  the  continued  fraction 
representation  of  ratios  of  Bessel  functions,  the  method  is  easy  to  use  and 
free  of  most  of  the  weaknesses  of  the  other  algorithms  used  to  generate  ratios 
of  Bessel  functions.2 

Let  us  review  the  method  of  computing  a  continued  fraction  by  defining  the 
simple  continued  fraction.  The  simple  continued  fraction  is  defined  for 
arbitrary  an  not  equal  to  zero  as 

F  =  a  +  I*° 

F  l  a 


rnr 


a'3  VI  ID 


a4  + 


(1) 


which  may  be  written  as 


(2) 


F  ,  a  +  2_  i_  J_ 

1  a2  +  a3  +  a4  + 


and  the  n'th  partial  convergent  of  the  continued  fraction  is 


‘Lentz,  W.  J.,  "Generating  Bessel  Functions  In  Mie  Scattering  Calculations 
Using  Continued  Fractions,”  Applied  Optics,  15:668.,  1976. 

‘Jaaskelalnen,  T.,  and  0.  Ruuskanen,  "Note  on  Lentz's  Algorithm,"  Applied 
Optics,  20:19,  1981. 
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The  Lentz  algorithm  Is  simply  a  way  of  confuting  the  Fn  without  sta^tlr-g  at 
some  indeterminate  last  an.  A  further  simplification  of  notation  is  needed: 

Fn  «  l*v  a2,  ....  an]  .  {4) 


Then  the  Lentz  algorithm  becomes 

Ca^]  [a2,  a^]  [aj,  a2*  •••  Ca^,  a^_^,  ...,  a^] 

Fn  U2J  La3,  a2J  ...  Lan’  an-l*  *  ’  a?J  ’  (5* 


Although  the  formula  looks  complicated,  the  successive  numerators  and 
denominators  are  generated  by  taking  the  reciprocal  of  the  previous  numerator 
or  denominator  and  adding  it  to  the  next  a^ 

When  a  given  numerator  is  equal  to  its  corresponding  denominator  (to  the 
number  of  digits  accuracy  that  is  desired),  the  continued  fraction  is  said  to 
have  converged. 


Original  Algorithm  Improvement 

The  original  algorithm  paper  included  an  improvement  technique  to  avoid 
inaccuracy  when  any  term  in  the  numerator  or  denominator  approached  or  became 
zero.  The  problem  step  is  bypassed  without  loss  of  accuracy  even  if  the  term 
is  exactly  zero.  Consider  the  following  case  in  which  1/a  is  equal  to  but  of 
opposite  sign  to  am.  When  am  is  added  to  1/a,  every  digit  that  is  the  same  in 

the  two  terms  will  cancel  and  cause  a  loss  of  accuracy  in  the  final  number  if 
no  steps  are  taken  to  bypass  the  problem.  For  example, 


Let 
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where  t  is  not  necessarily  1.0. 

The  product  of  the  two  terms  6y  =  t  can  be  correctly  calculated  by  equation 
(6).  This  avoids  the  problem  of  multiplying  zero  times  infinity  in  equation 
(5). 


t  is  the  product  of  the  next  two  terms.  z,  the  next  numerator  or  denominator, 
is  then 


z 


m+2 


Y  +  1  .  am+2  (Vl  6  +  U  +  8 


Vl  b"+T 


-  am+2  +  8/5  ' 


(7) 


The  algorithm  improvement  may  be  applied  to  the  denominator  by  replacing  ai  by 
a2  in  the  above  equations.  The  algorithm  improvement  may  be  applied 
independently  to  either  the  numerator  or  denominator  or  to  both. 

Numerical  Examples 

Consider  the  following  three  cases,  which  illustrate  the  implementation  of  the 
normal  Lentz  algorithm  and  the  Lentz  algorithm  with  improvement. 


Case  A-  F  =  /~Z  =  [1 ,2, 2, 2,2, . . .  ] 
na 

=  1*  3*  2.333*2.428571429*  2.411764706*  2.414634146 
r*’2. 5  *  2.4"*  2.4155666666  *  2.413793103 - 

Case  B:  Fnb  =  [1,1, -1,1,2]  =  0.5 
Case  C:  Fnc  =  [1 ,-l ,1 ,-l ,2]  =  0.5 


In  case  B  the  denominator  [-1,1]  is  exactly  zero,  regardless  of  the  computer 
word  length.  In  case  C  both  the  numerator  [-1,1]  and  the  denominator  [-1,1] 
are  zero,  causing  both  the  numerator  product  and  the  denominator  product  to  be 
zero.  For  case  B  one  has: 


.  [1]  [1,1]  [-1,1,1]  [1,-1, 1,1]  [2, 1,-1, 1,1] 
nb  - [11  [-1,1]  11,-1,11  [2,1, -1,1] - 


Applying  the  algorithm  improvement  one  has 


where  (1*0+1)  corresponds  to  equation  (6)  and  {2+0/1}  corresponds  to  equation 
(7).  The  continued  fraction  is  then  calculated  with  ease: 


1*2*  -0.5  *-1*1 

"-r+"r*"2 - 


0.5 


t 


In  like  manner,  the  improvement  may  be  applied  to  more  complicated  cases,  such 
as  case  C,  where  corrections  in  both  the  numerator  and  denominator  are 
requi red: 


.  Cl]  [-1,1]  Cl ,-1,1]  [-1,1, -1,1]  [2, -1,1, -1,1] 
nc  [-1]  II" -1]  r-U.'-lJ  L2, -171,-1] - 


In  this  case  we  have 


1  *  (1*0+1)  (-1+0/1)  (2+1/-1 ) 

r-n  T-PD+I)  C2+071) - 


0.5  . 


Notice  that  the  term  c  cannot  be  zero  unless  am+3  ■  0,  which  is  not  allowed  by 

definition.  If  any  a^  were  zero  In  the  derivation  of  a  representation  of  a 

function,  the  fraction  could  easily  be  transformed  to  the  form  of  equation  (2) 
by  formula  (5)  of  section  3.10.1  of  Abramowltz.’ 

Up  to  this  point  we  have  not  considered  the  Bessel  functions  for  which  the 
Lentz  algorithm  was  derived.  Like  Florida  orange  juice  for  breakfast,  the 
algorithm  is  not  just  for  Bessel  functions  anymore.  Rather,  it  is  a  general 
way  of  computing  the  continued  fraction  representation  of  any  function.  A 
good  source  of  continued  fraction  representations  of  functions  is  found  in 
Abramowltz  together  with  the  general  properties  of  the  functions. 

DISCUSSION  OF  THE  SIMPLIFICATION 

Simplified  Error  Removal 

In  some  cases  It  may  not  be  desirable  to  include  an  error  correction  when  e 
approaches  zero  but  Is  not  Identically  zero.  When  either  the  numerator  or 
denominator  Is  zero  to  the  accuracy  of  the  computer,  a  significant 
simplification  In  programming  is  possible.  This  simplification  reduces  the 


‘Abramowltz,  M.  and  I.  A.  Stegun,  ed.,  Handbook  of  Mathematical  Functions,  US 
Department  of  Commerce,  National  Bureau  of  Standards,  Washington,  DC,  20A02. 
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complexity  of  the  algorithm  and  will  often  speed  up  its  execution  with  little 
loss  of  accuracy. 

Consider  equations  (6)  and  (7)  when  8  is  exactly  zero: 

[am,. .. ,a23  =  B  =  0 

tam+l . al-*  =  Y  =  am+l  +  1*0/s 

tam+2 . al]  =  2  =  am+2  +  S/c  =  am+2  *  (8) 

And  the  product  of  the  two  terms  By  is 
Z  =  by  =  a^B  +  1  =  1. 

To  include  the  simplified  algorithm  improvement,  simply  do  not  take  the 
reciprocal  of  the  previous  numerator  and  denominator  before  adding  the  next 
am+2*  The  runmn9  product  is  multiplied  by  one,  so  that  it  does  not  change  if 

the  error  bypass  is  implemented.  This  simplification  was  first  published  in 
the  Hewlett  Packard  67  User's  Library  under  Spherical  Bessel  functions.*  The 
User's  Library  program  is  no  longer  available,  but  the  algorithm  in  appendix  A 
details  the  methods  used  in  that  program. 

SUMMARY 

The  original  algorithm  improvement  in  the  computation  of  continued  fractions 
may  be  applied  to  both  the  numerator  and  denominator  of  equation  (5).  A 
simplified  version  of  the  improvement  may  be  implemented  in  the  case  of  a 
numerator  or  denominator  being  exactly  zero.  A  FORTRAN  IV  program 
implementing  these  ideas  is  listed  for  ease  in  use. 


‘Lentz,  W.  J.,  Hewlett  Packard  67/97  User's  Library,  No  00642. 


APPENDIX  A 

FORTRAN  IV  PROGRAM  SIMPLE. FR 


Algorithm  Example 

The  implementation  of  the  simplified  form  of  the  algorithm  improvement  given 
in  equations  (6)  and  (7)  can  be  elusive.  The  best  method  to  avoid  pitfalls  is 
to  check  the  algorithm  against  well-defined  cases  such  as  Case  B  and  Case  C. 
The  FORTRAN  IV  program  Simple. FR,  listed  in  appendix  A,  correctly  calculates 

cases  A,  B,  and  C.  The  program  was  implemented  on  a  Data  General  Nova  3 

minicomputer  using  FORTRAN  V  revision  6.11.  The  nonstandard  accept  and  type 
statements  allow  input  and  output  in  arbitrary  format  and  may  be  replaced  by 
read  and  write  statements  in  standard  FORTRAN  IV.  The  variables  have  the 
following  functions: 

AN  =  the  input  terms  of  the  continued  fraction  from  equation  (4) 

NUM  =  the  current  partial  numerator  in  equation  (5) 

DEN  =  the  current  partial  denominator  in  equation  (5) 

PDT  =  the  n'th  convergent  of  the  continued  fraction  formed  1  the  n 
numerators  and  n-1  denominators  in  equation  (5) 

NFLAG  =  switch  to  avoid  taking  the  reciprocal  of  zero  In  th«  ,  rator 

DFLAG  =  switch  to  avoid  taking  the  reciprocal  of  zero  in  the  denominator 

ICOUNT  =  counter  to  keep  track  of  the  current  AN 


; 


SIKPLE.FR 


1 

2 

3 

4 


4 

10 
1  1 
12 

13 

14 

15 
Ur 
12 
13 
13 
20 
21 


25 

2  b 

7*0 

_  O 

29 

30 

31 

b2 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 


C  FROGRAM  S I r-lPLE  IS  DESIGNED  TO  HOVE  SEPARATE  L0OPS  P0P  NUM  AND  DEN 

1  A  SIMPLIFIED  TEST  IS  .  Vi 7  FOR  MJMERA7CR  MUD  DENOMINATOR  E0jA_  rr  ZEP9 

2  PROGRAM  PASSES  THE  TEST  OF  THE  FOLLOWING  SERIES 

•-  C  1.  1,-1.  1.21=0.5  AT  THE  uA3T  CONVERGENCE 

u  i  1.-1.  1. -1.23=0. 5  AT  IKE  LAST  CONVERGENCE 

C  THE  FIRST  TESTS  THE  DENOMINATOR  t  TRASS ,  AMD  THE  SECOND  THE  NL'lTE^pTQC 

1  BYPASS. 

L****#***:-;  ITT  ;  #.*  *v*W  *.  y «**■*■«  >  •■■**.  •  *.(  •*:» 

REAL  PD T, MUM. DEN, AH 

f.,.:l'**:*"'l-***  I  •;  "RROGf  m!  1  START  T  t  :tr 

I  TYPE  "ICOUNT-O,  INPUT  AH" 

ACCEPT  AN 

PDT=AH 
NO  1 1-  1  ./AH 
DEN  =  0 .  Q 
I l CUN  f-0 
HFLAG "0 
DEL  AG -i.; 

C****y*sm-:<  itlAIM  LOOP  POINT  AFTER  IN  ITIAL  IZATiON*****##**>toK*******#*#*******w** 

10  IFUIUM.RO. DEN. f4ND.NFLAG.NE.  1 .  AND .  DFLAG .  ME .  1 )  TYPE  "FRACTION  HAS  CONVERGED" 

!  COUNT--’  I COUNT >1 
TYPE  "PDT-  l!,PDT,  ICOUHT 
TYPE  "INPUT  AN" 

ACCEPT  AN 
IFiAN.EO.0)  GOTO  1 

C  STOP  IF  AN  IS  ZERO  UHICH  IS  NOT  ALLOWED  IN  THE  ALGORITHM 

C*****TEST  FOR  PREVIOUS  NUMERATOR  OR  DENOMINATOR  ZERO.  IF  SO  SHIP  STEP******** 

IF  INFLAG.NE. 1)  NUMUIUH+AN 
IF  (DFLAG. ME.  1)  DEN'-DEH-rAN 
NFLAP=NFI.AG+1 
DFLAG =DFLAC>1 

C*******’t:*.-!:rY  **;:OLLOW!NG  IS  NORMAL  NUMERATOR  LOOP  UNLESS  7ERn***************** 

IF (MUM. EO . 0)  GOTO  11 

PDT=PDT*NUU 

HUM- 1 . 0/NUI I 

NFL AG =0 

C*"!:***iK****i:**i:*!i:*FOLLOUING  IS  DENOMINATOR  LOOP'"******************************* 

II  IF(DEN.EO.O)  GOTO  10 
PDT  =PDT/1'EI! 

DEN  - 1 . 0/DEN 
DFLAG =0 
GOTO  10 
END 
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PO  Box  25007Denver  Federal  Center,  Bldg.  67 

Denver,  CO  80225 

Hugh  W.  Albers  (Executive  Secretary) 

CAO  Subcommittee  on  Atmos  Rsch 
National  Science  Foundation  Room  510 
Washington,  DC  2055 
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Or.  Eugene  W.  Bierly 
Director,  Division  of  Atmos  Sciences 
National  Scinece  Foundation 
1800  G  Street,  N.W, 

Washington,  DC  20550 

Commanding  Officer 
Naval  Research  Laboratory 
Code  2627 

Washington,  DC  20375 

Defense  Communications  Agency 
Technical  Library  Center 
Code  222 

Washington,  DC  20305 
Director 

Naval  Research  Laboratory 
Code  5530 

Washington,  X  20375 

Dr.  J.  M.  MacCallum 
Naval  Research  Laboratory 
Code  1409 

Washington,  X  20375 

HQDA  (DAEN-RDM/Dr.  de  Percin) 
Washington,  X  20314 

The  Library  of  Congress 
ATTN:  Exchange  A  Gift  Div 
Washington,  DC  20540 
2 

Mil  Asst  for  Atmos  Sci  Ofc  of 
the  Undersecretary  of  Defense 
for  Rsch  &  Engr/EALS  -  RM  3D129 
The  Pentagon 
Washington,  X  20301 

AFATL/DLODL 
Technical  Library 
Eglin  AFB,  FL  32542 

Naval  Training  Equipment  Center 
ATTN:  Technical  Information  Center 
Orlando,  FL  32813 

Technical  Library 

Chemical  Systems  Laboratory 

Aberdeen  Proving  Ground,  MD  21010 


US  Army  Materiel  Sy stans 
Analysis  Activity 
ATTN:  DRXSY-MP 
APG,  MD  21005 

Commander 

ERADCOM 

ATTN:  DRDEL-PA/ILS/-ED 
2800  Powder  Mill  Road 
Adelphi,  ND  20783 

Commander 

ERADCOM 

ATTN:  DRDEL-ST-T  (Dr.  B.  Zarwyn) 
2800  Powder  Mill  Road 
Adelphi,  MD  20783 
02 

Commander 

Harry  Diamond  Laboratories 
ATTN:  DELHD-CO 
2800  Powder  Mill  Road 
Adelphi,  MD  20783 

Chief 

Intel  Mat  Dev  A  Spt  Ofc 
ATTN:  DELEW-WL-I 
Bldg  4554 

Fort  George  G.  Mead,  MD  20755 

Acquisitions  Section,  IRDB-D823 
Library  A  Info  Svc  Div,  NOAA 
6009  Executive  Blvd. 

Rockville,  MD  20752 

Naval  Surface  Weapons  Center 
White  Oak  Library 
Silver  Spring,  MD  20910 

Air  Force  Geophysics  Laboratory 
ATTN:  LCC  (A.  S.  Carten,  Jr.) 
Hanscom  AFB,  MA  01731 

Air  Force  Geophysics  Laboratory 
ATTN:  LYD 

Hanscom  AFB,  MA  01731 

Meteorology  Division 
AFGL/LY 

Hanscom  AFB,  MA  D1731 
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The  Environmental  Research 
Institute  of  MI 
ATTN:  I R I A  Library 
PO  Box  8618 
Ann  Arbor,  MI  48107 

Mr.  will i am  A.  Main 
USDA  Forest  Service 
1407  S.  Harrison  Road 
East  Lansing,  MI  48823 

Dr.  A.  D.  Belmont 
Research  Divi sion 
PO  Box  1249 
Control  Data  Corp 
Minneapolis,  MN  55440 

Commander 

Naval  Oceanography  Command 
Bay  St.  Louis,  MS  39529 

Commanding  Officer 
US  Army  Armament  RAD  Command 
ATTN;  DRDAR-TSS  Bldg  59 
Dover,  NJ  07801 

Commander 

ERADCOM  Scientific  Advisor 

ATTN:  DRDEL-SA 

Fort  Monmouth,  NJ  07703 

Commander 

ERADCOM  Tech  Support  Activity 

ATTN:  DELSD-L 

Fort  Monmouth,  NJ  07703 

Commander 

HQ,  US  Army  Avionics  RAD  Ac  tv 

ATTN:  DAVAA-0 

Fort  Monmouth,  NJ  07703 

Commander 

USA  Elect  Warfare  Lab 
ATTN:  DELEW-DA  (File  Cy) 

Fort  Monmouth,  NJ  07703 

Commander 

US  Army  Electronics  RAD  Command 

ATTN:  DELCS-S 

Fort  Monmouth,  NJ  07703 


Commander 

US  Army  Satellite  Comm  Agency 

ATTN:  DRCPM-SC-3 

Fort  Monmouth,  NJ  07703 

Commander/Di  rector 
US  Army  Combat  Survl  A  Target 
Acquisition  Laboratory 
ATTN:  DELCS-0 
Fort  Monmouth,  NJ  07703 

Oi rector 

Night  Vision  A  Electro-Optics  Laboratory 
ATTN:  DELNV-L  (Dr.  R.  Buser) 

Fort  Bel  voir,  VA  22060 

Project  Manager 
FIREFINDER/REMBASS 
ATTN:  DRCPM-FFR-TM 
Fort  Monmouth,  NJ  07703 

6585  TG/WE 

Holloman  AFB,  NM  88330 

AFWL/Technical  Library  (SUL) 

Kirtland  AFB,  NM  87117 

AFWl/WE 

Kirtland,  AFB,  NM  87117 
TRASANA 

ATTN:  ATAA-SL  (D.  Anguiano) 

WSMR,  NM  88002 

Commander 

US  Army  White  Sands  Missile  Range 
ATTN:  STEWS-PT-AL 

White  Sands  Missile  Range,  W  88002 

Rome  Air  Development  Center 
ATTN:  Documents  Library 
TSLD  (Bette  Smith) 

Griffiss  AFB,  NY  13441 

Environmental  Protection  Agency 
Meteorology  Laboratory,  MD  80 
Rsch  Triangle  Park,  NC  27711 


US  Army  Research  Office 
ATTN:  DRXRO-PP 
PO  Box  12211 

Rsch  Triangle  Park,  NC  27709 
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Commandant 

US  Army  Field  Artillery  School 
ATTN:  ATSF 'CO-MS  (Mr.  Farmer) 

Fort  Sill,  OK  73503 

Commandant 

US  Army  Field  Artillery  School 
ATTN:  ATSF-CF-R 
Fort  Sill,  OK  73503 

Commandant 

US  Army  Field  Artillery  School 
ATTN:  Morris  Swett  Library 
Fort  Sill,  OK  73503 

Commander 

US  Army  Dugway  Proving  Ground 
ATTN:  STFDP-MT  -DA-M 

(Mr.  Paul  Carlson) 

Ouqway,  UT  84022 

Commander 

US  Army  Ouqway  Proving  Ground 
ATTN:  MT-DA-L 
Dugway,  UT  84022 

US  Army  Ougway  Proving  Ground 
ATTN:  STEDP-MT-DA-T 

(Or.  W.  A.  Peterson) 

Dugway,  UT  84022 

Inge  Dirmhirn,  Professor 
Utah  State  University,  UMC  48 
Logan,  UT  84322 

Defense  Technical  Information  Center 
ATTN:  DTIC-DDA-2 
Cameron  Station,  Bldg.  5 
Alexandria,  VA  22314 
12 

Commanding  Officer 

US  Army  Foreign  Sci  4  Tech  Cen 

ATTN:  DRXST-IS1 

220  7th  Street,  NE 

Charlottesville,  VA  22901 

Naval  Surface  weapons  Center 
Code  G65 

Dahlgren,  VA  22448 


Commander 

US  Army  Night  Vision 
S  Electro-Optics  Lab 
ATTN:  OELNV-D 
Fort  Bel  voir,  VA  22060 

Commander 
USATRADOC 
ATTN:  ATCO-FA 
Fort  Monroe,  VA  23651 

Commander 
USATRADOC 
ATTN:  ATCD-IR 
Fort  Monroe,  VA  23651 

Dept  of  the  Air  Force 
5WW/DN 

Langley  AF8,  VA  23665 

US  Army  Nuclear  S  Cml  Agency 
ATTN:  MONA-WE 
Springfield,  VA  22150 

Di rector 

US  Army  Signals  Warfare  Lab 
ATTN:  DELSW-OS  (Dr.  Burkhardt) 
Vint  Mill  Farms  Station 
Warrenton,  VA  22186 

Commander 

US  Army  Cold  Regions  Test  Cen 
ATTN:  STECR-OP-PM 
APO  Seattle,  WA  98733 
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